/* Copyright 2007-2008 dnAnalytics Project.
 *
 * Contributors to this file:
 * Patrick van der Velde
 *
 * This file is part of dnAnalytics.  dnAnalytics is licensed under the 
 * Microsoft Public License. See License.txt for a complete copy of the
 * license.
 */

using System;
using System.Diagnostics;
using dnAnalytics.LinearAlgebra.Solvers.Iterative;
using dnAnalytics.LinearAlgebra;
using NUnit.Framework;
using dnAnalytics.Math;
using dnAnalytics.LinearAlgebra.Solvers;

namespace dnAnalytics.Tests.LinearAlgebra.Solvers.Iterative
{
    [TestFixture]
    [Category("Managed")]
    public class MlkBiCgStabTest
    {
        private const double s_ConvergenceBoundary = 1e-10;
        private const int s_MaximumIterations = 1000;

        [Test]
        [ExpectedException(typeof(MatrixNotSquareException))]
        public void SolveWideMatrix()
        {
            SparseMatrix matrix = new SparseMatrix(2, 3);
            Vector input = new DenseVector(2);

            MlkBiCgStab solver = new MlkBiCgStab();

            Vector output = solver.Solve(matrix, input);
        }

        [Test]
        [ExpectedException(typeof(MatrixNotSquareException))]
        public void SolveLongMatrix()
        {
            SparseMatrix matrix = new SparseMatrix(3, 2);
            Vector input = new DenseVector(3);

            MlkBiCgStab solver = new MlkBiCgStab();

            Vector output = solver.Solve(matrix, input);
        }

        [Test]
        public void SolveUnitMatrixAndBackMultiply()
        {
            // Create the identity matrix
            Matrix matrix = new SparseMatrix(100);
            for (int i = 0; i < matrix.Rows; i++)
            {
                matrix[i, i] = 1.0;
            }

            // Create the y vector
            Vector y = new DenseVector(matrix.Rows, 1);

            // Create an iteration monitor which will keep track of iterative convergence
            Iterator monitor = new Iterator(new IIterationStopCriterium[]
                {
                    new IterationCountStopCriterium(s_MaximumIterations),
                    new ResidualStopCriterium(s_ConvergenceBoundary),
                    new DivergenceStopCriterium(),
                    new FailureStopCriterium()
                });
            MlkBiCgStab solver = new MlkBiCgStab(monitor);

            // Solve equation Ax = y
            Vector x = solver.Solve(matrix, y);

            // Now compare the results
            Assert.IsNotNull(x, "#02");
            Assert.AreEqual(y.Count, x.Count, "#03");

            // Back multiply the vector
            Vector z = matrix.Multiply(x);

            // Check that the solution converged
            Assert.IsTrue(monitor.Status is CalculationConverged, "#04");

            // Now compare the vectors
            for (int i = 0; i < y.Count; i++)
            {
                Assert.IsTrue(Precision.IsSmallerWithTolerance(y[i] - z[i], s_ConvergenceBoundary, 1), "#05-" + i.ToString());
            }
        }

        [Test]
        public void SolveScaledUnitMatrixAndBackMultiply()
        {
            // Create the identity matrix
            Matrix matrix = new SparseMatrix(100);
            for (int i = 0; i < matrix.Rows; i++)
            {
                matrix[i, i] = 1.0;
            }
            // Scale it with a funny number
            matrix.Multiply(System.Math.PI);

            // Create the y vector
            Vector y = new DenseVector(matrix.Rows, 1);

            // Create an iteration monitor which will keep track of iterative convergence
            Iterator monitor = new Iterator(new IIterationStopCriterium[]
                {
                    new IterationCountStopCriterium(s_MaximumIterations),
                    new ResidualStopCriterium(s_ConvergenceBoundary),
                    new DivergenceStopCriterium(),
                    new FailureStopCriterium()
                });
            MlkBiCgStab solver = new MlkBiCgStab(monitor);

            // Solve equation Ax = y
            Vector x = solver.Solve(matrix, y);

            // Now compare the results
            Assert.IsNotNull(x, "#02");
            Assert.AreEqual(y.Count, x.Count, "#03");

            // Back multiply the vector
            Vector z = matrix.Multiply(x);

            // Check that the solution converged
            Assert.IsTrue(monitor.Status is CalculationConverged, "#04");

            // Now compare the vectors
            for (int i = 0; i < y.Count; i++)
            {
                Assert.IsTrue(Precision.IsSmallerWithTolerance(y[i] - z[i], s_ConvergenceBoundary, 1), "#05-" + i.ToString());
            }
        }

        [Test]
        public void SolvePoissonMatrixAndBackMultiply()
        {
            // Create the matrix
            SparseMatrix matrix = new SparseMatrix(100);
            // Assemble the matrix. We assume we're solving the Poisson equation
            // on a rectangular 10 x 10 grid
            int gridSize = 10;
            // The pattern is:
            // 0 .... 0 -1 0 0 0 0 0 0 0 0 -1 4 -1 0 0 0 0 0 0 0 0 -1 0 0 ... 0
            for (int i = 0; i < matrix.Rows; i++)
            {
                // Insert the first set of -1's
                if (i > (gridSize - 1))
                {
                    matrix[i, i - gridSize] = -1;
                }

                // Insert the second set of -1's
                if (i > 0)
                {
                    matrix[i, i - 1] = -1;
                }

                // Insert the centerline values
                matrix[i, i] = 4;

                // Insert the first trailing set of -1's
                if (i < matrix.Rows - 1)
                {
                    matrix[i, i + 1] = -1;
                }

                // Insert the second trailing set of -1's
                if (i < matrix.Rows - gridSize)
                {
                    matrix[i, i + gridSize] = -1;
                }
            }

            // Create the y vector
            Vector y = new DenseVector(matrix.Rows, 1);

            // Create an iteration monitor which will keep track of iterative convergence
            Iterator monitor = new Iterator(new IIterationStopCriterium[]
                {
                    new IterationCountStopCriterium(s_MaximumIterations),
                    new ResidualStopCriterium(s_ConvergenceBoundary),
                    new DivergenceStopCriterium(),
                    new FailureStopCriterium()
                });
            MlkBiCgStab solver = new MlkBiCgStab(monitor);

            // Solve equation Ax = y
            Vector x = solver.Solve(matrix, y);

            // Now compare the results
            Assert.IsNotNull(x, "#02");
            Assert.AreEqual(y.Count, x.Count, "#03");

            // Back multiply the vector
            Vector z = matrix.Multiply(x);

            // Check that the solution converged
            Assert.IsTrue(monitor.Status is CalculationConverged, "#04");

            // Now compare the vectors
            for (int i = 0; i < y.Count; i++)
            {
                Assert.IsTrue(Precision.IsSmallerWithTolerance(System.Math.Abs(y[i] - z[i]), s_ConvergenceBoundary, 1), "#05-" + i.ToString());
            }
        }
    }
}
